########################
## MEDIATION ANALYSIS ##
## 29 OCT 2019        ##
########################

## LIBRARY----------------------------------------------------------------------
library(readstata13) #read dta
library(mediation) #mediation analysis

## READ IN ---------------------------------------------------------------------

sparcs_df <- read.dta13('sparcs.dta')

#only want respondents who were present in both rounds
sparcs_df <- sparcs_df[which(sparcs_df$both_rounds == 1), ]

#set seed for reproducibility
set.seed(1017)

#subset to only needed variables, standardize
rel_vars <- c('putin_approval', 'corrupt_high', 'corrupt_low', 'rusecon_fut', 
              'changes_fam', 'state_tv_cats', 'discuss_ind', 'emot_pos_ind',
              'live', 'org_empl', 'fin_stat', 'educ', 'cohort', 'famecon', 
              'PanelResp')

stand_df <- sparcs_df[ , rel_vars]

stand_df <- apply(stand_df, 2, scale)
stand_df <- as.data.frame(stand_df)


## FUNCTION--------------------------------------------------------

mmed.function <- function(treat, med, outcome, df, out_dir){

  stopifnot(is.character(treat) & is.character(med) & is.character(outcome))
  
  #need to subset to observations not missing on treatment
  nas <- which(is.na(df[ , outcome]))
  
  if(length(nas) > 0) {
  df <- df[-nas, ]
  }
  #if no missing - can use whole df
  
  #mediator and outcome models
  med.mod <- lm(eval(parse(text = med)) ~ eval(parse(text = treat)) + live + org_empl + fin_stat + educ + cohort + famecon + 
                      factor(df$PanelResp), 
                    data = df) #mediator model
  
  out.mod <- lm(eval(parse(text = outcome))~ eval(parse(text = treat)) + eval(parse(text = med))
                  + live + org_empl + fin_stat 
                  + educ + cohort + famecon + factor(df$PanelResp),
                  data = df) #outcome model
  
  fe.res <- mediate(model.m = med.mod, 
                         model.y = out.mod, 
                         treat =  "eval(parse(text = treat))",
                         mediator = "eval(parse(text = med))",
                         sims = 1000,
  )
  

  #write out
  sink(file = paste0(out_dir, treat, outcome, ".txt"))
  print(summary(fe.res))
  sink()
  
  closeAllConnections()
  
  print(paste(treat, outcome, "done"))
  
}

# TABLE 2-----------------------------------------------------------------------

outcomes <- c('putin_approval', 'corrupt_high', 'corrupt_low', 'rusecon_fut', 
              'changes_fam')

sapply(outcomes, mmed.function, treat = 'state_tv_cats', med = 'emot_pos_ind',
       out_dir = './table_2/', df = stand_df)

# APP---------------------------------------------------------------------------

#discuss_ind
sapply(outcomes, mmed.function, treat = 'discuss_ind', med = 'emot_pos_ind',
       out_dir = './appendix/', df = sparcs_df)

#polnews
sapply(outcomes, mmed.function, treat = 'polnews', med = 'emot_pos_ind',
       out_dir = './appendix/', df = sparcs_df)


